﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using MPI;

namespace Dwarf.Managed
{
    /// <summary>
    /// Solve the Structured Grid problem with MPI parallelization  
    /// </summary>
    /// <remarks> 
    /// For the demonstration purposes of the project, we have chosen solution of Laplace's equation 
    /// with the Dirichlet’s conditions on the square. To obtain solution, iterated Jacobi method is applied. 
    /// It consists of steps like:
    /// newgrid[i,j] = (grid[i-1,j] + grid[i+1,j] + grid[i,j - 1] + grid[I,j + 1]) / 4
    /// </remarks>
    class Solver
    {
        /// <summary>
        /// The number of columns.
        /// </summary>
        public int xmax;

        /// <summary>
        /// The number of rows.
        /// </summary>
        public int ymax;

        /// <summary>
        /// Epsilon.
        /// </summary>
        public double epsilon;

        /// <summary>
        /// Initial matrix.
        /// </summary>
        public double[][] grid;

        /// <summary>
        /// Iterations taken
        /// </summary>
        public int iterations;

        /// <summary>
        /// Temp matrix.
        /// </summary>
        public double[][] tempGrid;

        /// <summary>
        /// Empty constructor.
        /// </summary>
        public Solver()
        {
        }

        /// <summary>
        /// MPI based method.
        /// </summary>
        public void Solve()
        {
            int countRows = Communicator.world.Rank == Settings.ROOT_PROCESS ?      //Number of rows.
                (ymax - 2) / Communicator.world.Size + 2 + (ymax - 2) % Communicator.world.Size :
                (ymax - 2) / Communicator.world.Size + 2;

            bool converged = false;                                                 //Convergence of matrix
            iterations = 0;

            while (!converged)                                                      //Main loop.
            {
                iterations++;
                converged = true;

                if (Communicator.world.Rank == 0)
                {
                    //Send values from 0 process to 1
                    Communicator.world.Send(grid[countRows - 2], Communicator.world.Rank + 1, Communicator.world.Rank + 1);
                    //Receive values from 1 process by 0
                    Communicator.world.Receive(Communicator.world.Rank + 1, Communicator.anyTag, ref grid[countRows - 1]);
                }
                else if (Communicator.world.Rank == Communicator.world.Size - 1)
                {
                    //Receive values by last process
                    Communicator.world.Receive(Communicator.world.Rank - 1, Communicator.anyTag,ref grid[0]);
                    //Send values from last process
                    Communicator.world.Send(grid[1], Communicator.world.Rank - 1, Communicator.world.Rank - 1);
                }
                else
                {
                    //Send & receive values for other rows
                    Communicator.world.Receive(Communicator.world.Rank - 1, Communicator.anyTag,ref grid[0]);
                    Communicator.world.Send(grid[1], Communicator.world.Rank - 1, Communicator.world.Rank - 1);
                    Communicator.world.Send(grid[countRows - 2], Communicator.world.Rank + 1, Communicator.world.Rank + 1);
                    Communicator.world.Receive(Communicator.world.Rank + 1, Communicator.anyTag,ref grid[countRows - 1]);
                }

                for (int i = 1; i < countRows - 1; i++)                                                                          // Loop for number of rows in matrix.
                {
                    //Since C# compiler/runtime is not hoisting 
                    double[] above = grid[i - 1];
                    double[] row = grid[i];
                    double[] below = grid[i + 1];
                    double[] tempRow = tempGrid[i];
                    for (int j = 1; j < xmax - 1; j++)                                                                      // Loop for number of columns in matrix.
                    {
                        //Perform step of computation on current node
                        double cell = tempRow[j] = (above[j] + below[j] + row[j - 1] + row[j + 1]) / 4.0;

                        converged = (converged && (Math.Abs(cell - row[j]) < epsilon));         // Check test of convergence for current node                     }
                    }
                }

                //Using MPI_Allreduce() to distribute convergence
                converged = Communicator.world.Allreduce(converged, MPI.Operation<bool>.LogicalAnd);

                //swap newly computed tempGrid and grid to get ready for next interation
                double[][] swap = grid;
                grid = tempGrid;
                tempGrid = swap;

                Communicator.world.Barrier();
            }

            //Collect the resulting node’s values on root process
            if (Communicator.world.Rank == Settings.ROOT_PROCESS)
            {
                int count = (ymax - 2) / Communicator.world.Size + 2;
                for (int i = 1; i < Communicator.world.Size; i++)
                {
                    for (int j = 1; j < count - 1; j++)
                        Communicator.world.Receive(
                            i, 
                            i + j, 
                            ref grid[(ymax - 2) % Communicator.world.Size + ((ymax - 2) / Communicator.world.Size) * i + j]);
                }
            }
            else
            {
                for (int i = 1; i < countRows - 1; i++)
                    Communicator.world.Send(
                        grid[i],
                        Settings.ROOT_PROCESS,
                        Communicator.world.Rank + i);
            }

            //grid = (odd ? grid : tempGrid);            
        }
    }
}
